MODULE menu_prim
  !-----------------------------------------------------------
  ! KRUUNUNHAKA basis set tool kit - primitive optimization menu
  !-----------------------------------------------------------
  ! See CHANGELOG
  !-----------------------------------------------------------
  ! Copyright (C) 2003-2008 Pekka Manninen, 2010 Jussi Lehtola
  !
  ! This program is distributed under the terms of the GNU General
  ! Public License as published by the Free Software Foundation;
  ! either version 2 of the License, or (at your option) any later version.
  !
  ! This program is distributed in the hope that it will be useful,
  ! but WITHOUT ANY WARRANTY; without even the implied warranty of
  ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. 
  !
  ! Please note that this program as well as its 
  ! author must be properly cited whenever the program or some parts 
  ! originated on it are employed. 
  !
  !-----------------------------------------------------------
  USE definitions
  USE textutil
  USE menu
  USE resize
  USE optimize_eventempered
  USE optimize_primitives
  USE sorting
  IMPLICIT NONE
  
CONTAINS
 !-----------------------------------------------------------  
  SUBROUTINE prim_menu(X_l,C_l,nbas)
    IMPLICIT NONE
    REAL(KIND=prec), DIMENSION(:,:), ALLOCATABLE, INTENT(INOUT) :: X_l
    REAL(KIND=prec), DIMENSION(:,:,:), ALLOCATABLE, INTENT(INOUT) :: C_l
    INTEGER, DIMENSION(MAXL,2), INTENT(INOUT) :: nbas
    CHARACTER(LEN=TERMWIDTH) :: basmenu_tit
    CHARACTER(LEN=TERMWIDTH), DIMENSION(12) :: basmenu
    CHARACTER(LEN=12) :: bassel
    INTEGER :: sel
    ! Current values
    CHARACTER(LEN=TERMWIDTH), DIMENSION(4) :: currentvals
    ! Angular momentum to optimize for
    INTEGER, SAVE :: l=1
    ! Exponent range in logarithm scale!
    REAL(KIND=prec), SAVE :: low=-1.0_prec
    REAL(KIND=prec), SAVE :: high=1.0_prec
    ! Number of primitives
    INTEGER, SAVE :: nprim=5
    ! Measure of goodness to use
    INTEGER, SAVE :: measure=2
    ! Use STO to probe?
    LOGICAL, SAVE :: usesto=.FALSE.
    ! Loop index
    INTEGER :: ix
    LOGICAL :: finished

    ! Wanted tolerance
    REAL(KIND=prec), SAVE :: tol
    ! Current result
    REAL(KIND=prec) :: curres
    ! Maximum number of primitives
    INTEGER, PARAMETER :: maxprim=150
    ! Time variables
    INTEGER :: time1, time2, cpucount

    basmenu_tit="Primitive set optimization menu"
    basmenu(1)='Set angular momentum type'
    basmenu(2)='Set range to optimize in'
    basmenu(3)='Change measure of goodness of fit'
    basmenu(4)='Change type of probe'
    basmenu(5)='Show current primitives'
    basmenu(6)='Create (unoptimized) even-tempered set'
    basmenu(7)='Create optimized even-tempered set'
    basmenu(8)='Create optimized well-tempered set'
    basmenu(9)='Create fully optimized set by gradient walk'
    basmenu(10)='Create fully optimized set w/ evolution algorithm'
    basmenu(11)='Optimize wrt number of primitives to reach a given tolerance'
    basmenu(12)='Return to main menu'
    bassel='123456789esr'
    
    finished=.FALSE.

    DO WHILE(.not. finished)
       ! Form current values
       WRITE(currentvals(1),'(A,A,A,I3,A)') 'Working on the ',l_name(l),&
            ' shell currently containing ', nbas(l,1),' primitives.'
       WRITE(currentvals(2),'(A,F6.2,A,F6.2,A)') &
            'Trying to optimize completeness from ', low, ' to ', high, '.'
       IF(usesto) THEN
          WRITE(currentvals(3),'(A)') &
               'Using a Slater Type Orbital as the probe.'
       ELSE
          WRITE(currentvals(3),'(A)') &
               'Using a Gaussian Orbital as the probe.'
       END IF

       ! Which measure of goodness is used?
       SELECT CASE(measure)
       CASE(1)
          WRITE(currentvals(4),'(A)') 'Trying to minimize maximum deviance.'
       CASE(2)
          WRITE(currentvals(4),'(A)') &
               'Trying to minimize rms deviation from unity.'
       CASE(3)
          WRITE(currentvals(4),'(A)') &
               'Trying to maximize integral over profile.'
       END SELECT

       ! Do selection
       sel=make_cvals_selection(basmenu_tit,currentvals,basmenu,bassel,'*')
       SELECT CASE(sel)
       CASE(1)
          l=-1
          DO WHILE(l < 0 .OR. l > MAXL)
             l=read_int('Give value of angular momentum (S = 0, P = 1, etc)')
             IF (l < 0 .OR. l > MAXL) THEN
                WRITE (*,'(A)') 'Invalid angular momentum.'
             END IF
          END DO
          ! Convert to Fortran indexing
          l=l+1
       CASE(2)
          ! Set limits
          low=1.0_prec
          high=-1.0_prec
          DO WHILE(low>high)
             CALL read_limits('Give limits of scanning range in &
             &10-base logarithm scale:', low, high)

             IF(low>high) THEN
                WRITE (*,'(A)') 'Invalid scanning range.'
             END IF
          END DO
       CASE(3)
          ! Change measure of goodness
          measure=0
          DO WHILE(measure <1 .OR. measure>3)
             measure=read_int('1. Minimize deviance, 2. Minimize RMS deviance, 3. Maximize integral?')
             IF(measure <1 .OR. measure>3) THEN
                WRITE (*,'(A)') 'Invalid measure.'
             END IF
          END DO
       CASE(4)
          ! Change probe
          usesto=.NOT. usesto
       CASE(5)
          ! Show current primitives
          IF(nbas(l,1)>0) THEN
             WRITE (*,'(A,A,A)') 'Current primitives for the ',&
                  l_name(l),' shell are:'
             DO ix=1,nbas(l,1)
                WRITE (*,'(E20.8)') X_l(l,ix)
             END DO
          ELSE
             WRITE(*,'(A,A,A)') 'There are currently no primitives on the ',&
                  l_name(l), ' shell.'
          END IF

          WRITE (*,*) ''
          
       CASE(6)
          ! Just create even tempered set

          ! Check if shell already has something
          IF(clearok(nbas,l)) THEN
             ! Get wanted number of primitives
             CALL get_nprim(nprim)
             ! Check sufficient memory allocation
             CALL resize_Xlcl(nbas,X_l,C_l,nprim)
             ! Initialize exponents
             X_l(l,:)=0.0_prec
             nbas(l,1)=nprim
             CALL first_set(low, high, X_l(l,1:nbas(l,1)))
             ! Make dummy contraction matrix
             C_l(l,:,:)=0.0_prec
             DO ix=1,nprim
                C_l(l,ix,ix)=1.0_prec
             END DO
             nbas(l,2)=nprim
          END IF
       CASE(7)
          ! Perform optimization of even-tempered set

          IF(clearok(nbas,l)) THEN 
             ! Get wanted number of primitives
             CALL get_nprim(nprim)

             ! Check sufficient memory allocation
             CALL resize_Xlcl(nbas,X_l,C_l,nprim)
             
             ! Initialize exponents
             nbas(l,1)=nprim
             CALL first_set(10**low, 10**high, X_l(l,1:nbas(l,1)))
             
             ! Make dummy contraction matrix
             C_l(l,:,:)=0.0_prec
             DO ix=1,nprim
                C_l(l,ix,ix)=1.0_prec
             END DO
             nbas(l,2)=nprim
             
             ! Run optimization
             curres=opt_eventempered(l, 10**low, 10**high, X_l(l,1:nbas(l,1)),&
                  C_l(l,1:nbas(l,1),1:nbas(l,2)), usesto, measure, .TRUE.)

          END IF
       CASE(8)
          ! Perform optimization of well-tempered set

          IF(clearok(nbas,l)) THEN 
             ! Get wanted number of primitives
             CALL get_nprim(nprim)
             ! Check sufficient memory allocation
             CALL resize_Xlcl(nbas,X_l,C_l,nprim)
             ! Initialize exponents
             nbas(l,1)=nprim
             CALL first_set(10**low, 10**high, X_l(l,1:nbas(l,1)))
             
             ! Make dummy contraction matrix
             C_l(l,:,:)=0.0_prec
             DO ix=1,nprim
                C_l(l,ix,ix)=1.0_prec
             END DO
             nbas(l,2)=nprim
             
             ! Run optimization
             curres=opt_welltempered(l, 10**low, 10**high, X_l(l,1:nbas(l,1)),&
                  C_l(l,1:nbas(l,1),1:nbas(l,2)), usesto, measure, .TRUE.)
          END IF
       CASE(9)
          ! Perform optimization

          IF(clearok(nbas,l)) THEN 
             ! Get wanted number of primitives
             CALL get_nprim(nprim)
             ! Check sufficient memory allocation
             CALL resize_Xlcl(nbas,X_l,C_l,nprim)
             
             ! Initialize exponents
             nbas(l,1)=nprim
             CALL first_set(10**low, 10**high, X_l(l,1:nbas(l,1)))
             
             ! Make dummy contraction matrix
             C_l(l,:,:)=0.0_prec
             DO ix=1,nprim
                C_l(l,ix,ix)=1.0_prec
             END DO
             nbas(l,2)=nprim

             WRITE (*,'(A)',ADVANCE='no') 'Initializing with an optimized &
                  &even tempered set ...'
             curres=opt_eventempered(l, 10**low, 10**high, X_l(l,1:nbas(l,1)),&
                  C_l(l,1:nbas(l,1),1:nbas(l,2)), usesto, measure, .FALSE.)
             WRITE (*,'(A)') ' done.'

             WRITE (*,'(A)') 'Proceeding with full optimization.'
             curres=opt_prim_bf(l, 10**low, 10**high, X_l(l,1:nbas(l,1)), &
                  C_l(l,1:nbas(l,1),1:nbas(l,2)), usesto, measure,.TRUE.)
          END IF

       CASE(10)
          ! Perform optimization
          
          IF(clearok(nbas,l)) THEN 
             ! Get wanted number of primitives
             CALL get_nprim(nprim)
             ! Check sufficient memory allocation
             CALL resize_Xlcl(nbas,X_l,C_l,nprim)
             
             ! Initialize exponents
             nbas(l,1)=nprim
             CALL first_set(10**low, 10**high, X_l(l,1:nbas(l,1)))
             
             ! Make dummy contraction matrix
             C_l(l,:,:)=0.0_prec
             DO ix=1,nprim
                C_l(l,ix,ix)=1.0_prec
             END DO
             nbas(l,2)=nprim

             ! Run optimization
             CALL opt_prim(l, 10**low, 10**high, X_l(l,1:nbas(l,1)), &
                  C_l(l,1:nbas(l,1),1:nbas(l,2)), usesto, measure)
          END IF
          
       CASE(11)
          ! Determine number of functions necessary to reach a given tolerance
          ! (also get exponents)

          ! Read in wanted tolerance
          tol=-1.0_prec
          DO WHILE(tol<0.0_prec)
             tol=read_double('Give wanted tolerance:')
             IF(tol<0) THEN
                WRITE (*,'(A)') 'Invalid tolerance.'
             END IF
          END DO

          IF(clearok(nbas,l)) THEN
             ! Loop over amount of primitives
             WRITE (*,'(A25,A15,A15)') 'N','Measure','Time (s)'
             DO nprim=1,maxprim

                ! Get time
                CALL system_clock(time1,cpucount)
                
                ! Check sufficient memory allocation
                CALL resize_Xlcl(nbas,X_l,C_l,nprim)
                
                ! Initialize exponents
                nbas(l,1)=nprim
                CALL first_set(10**low, 10**high, X_l(l,1:nprim))
                
                ! Make dummy contraction matrix
                C_l(l,:,:)=0.0_prec
                DO ix=1,nprim
                   C_l(l,ix,ix)=1.0_prec
                END DO
                nbas(l,2)=nprim
                
                ! Initialize with even tempered set
                curres=opt_eventempered(l, 10**low, 10**high, &
                     X_l(l,1:nbas(l,1)), C_l(l,1:nbas(l,1),1:nbas(l,2)), &
                     usesto, measure, .FALSE.)
                ! and optimize it fully
                curres=opt_prim_bf(l, 10**low, 10**high, X_l(l,1:nbas(l,1)),&
                     C_l(l,1:nbas(l,1),1:nbas(l,2)), usesto, measure, .FALSE.)

                ! Get time
                CALL system_clock(time2,cpucount)

                ! Print out status
                WRITE(*,'(I25,E15.4,F15.2)') nprim, curres, &
                     (time2-time1)*1.0/cpucount

                ! Have we reached the tolerance already?
                IF(curres<tol) EXIT
             END DO
          END IF
       CASE(12)
          ! Return to main menu
          finished=.TRUE.
       END SELECT
    END DO
  END SUBROUTINE prim_menu
  !-----------------------------------------------------------
  FUNCTION clearok(nbas, l)
    ! Is it OK to discard existing primitives?
    IMPLICIT NONE
    INTEGER, INTENT(IN), DIMENSION(MAXL,2) :: nbas
    INTEGER, INTENT(IN) :: l
    ! Helper array
    CHARACTER(LEN=TERMWIDTH) :: tmpstr
    LOGICAL :: clearok

    clearok=.TRUE.
    
    IF(nbas(l,1)>0) THEN
       WRITE (tmpstr,'(A,A,A,I3,A)') 'The ', l_name(l), &
            ' shell already has ', nbas(l,1), ' primitives.'
       IF(make_yn(tmpstr)) THEN
          ! OK to proceed
          clearok=.TRUE.
       ELSE
          ! Not ok.
          clearok=.FALSE.
       END IF
    END IF
  END FUNCTION clearok
  !-----------------------------------------------------------
  SUBROUTINE get_nprim(nprim)
    ! Get wanted number of primitives
    IMPLICIT NONE
    INTEGER, INTENT(OUT) :: nprim
    
    nprim=-1
    DO WHILE(nprim <= 0)
       nprim=read_int('Give wanted amount of primitives:')
       IF(nprim <=0) THEN
          WRITE (*,'(A)') 'Invalid number of primitives.'
       END IF
    END DO
  END SUBROUTINE get_nprim
  !-----------------------------------------------------------
  SUBROUTINE first_set(low, up, init_set)
    ! Initialize exponents as even tempered set
    IMPLICIT NONE
    REAL(KIND=prec), DIMENSION(:), INTENT(OUT) :: init_set
    REAL(KIND=prec), INTENT(IN) :: low, up
    REAL(KIND=prec) :: alpha, beta

    INTEGER :: n
    n=SIZE(init_set)

    ! In an even tempered set the exponents are distributed according to
    ! zeta(i) = alpha * beta^(i-1)

    beta = EXP(LOG(up/low)/(REAL(n,prec)))
    alpha=low*SQRT(beta)

    CALL form_eventempered(alpha, beta, init_set)
  END SUBROUTINE first_set
  !----------------------------------------------------------- 
END MODULE menu_prim
